!! Copyright (C) 2009,2010,2011,2012  Marco Restelli
!!
!! This file is part of:
!!   FEMilaro -- Finite Element Method toolkit
!!
!! FEMilaro is free software; you can redistribute it and/or modify it
!! under the terms of the GNU General Public License as published by
!! the Free Software Foundation; either version 3 of the License, or
!! (at your option) any later version.
!!
!! FEMilaro is distributed in the hope that it will be useful, but
!! WITHOUT ANY WARRANTY; without even the implied warranty of
!! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
!! General Public License for more details.
!!
!! You should have received a copy of the GNU General Public License
!! along with FEMilaro; If not, see <http://www.gnu.org/licenses/>.
!!
!! author: Marco Restelli                   <marco.restelli@gmail.com>

!>\brief
!!
!! Laundau damping problem.
!!
!! \n
!!
!! For references, consider the following papers: 
!! <a
!! href="http://dx.doi.org/10.1016/S0010-4655(99)00247-7">[Nakamura,
!! Yabe, 1999]</a>, <a
!! href="http://dx.doi.org/10.1006/jcph.2001.6818">[Filbet,
!! Sonnendr&uuml;cker, Bertrand, 2001]</a>, <a
!! href="http://dx.doi.org/10.1016/S0010-4655(02)00694-X">[Filbet,
!! Sonnendr&uuml;cker, 2003]</a>,
!! <a
!! href="http://dx.doi.org/10.1051/proc/201343003">[Mehrenberger,
!! Steiner, Marradi, Crouseilles, Sonnendrücker, Afeyan, 2013]</a>.
!<----------------------------------------------------------------------
module mod_landau_damping_test

!-----------------------------------------------------------------------

 use mod_messages, only: &
   mod_messages_initialized, &
   error,   &
   warning, &
   info

 use mod_kinds, only: &
   mod_kinds_initialized, &
   wp

 use mod_fu_manager, only: &
   mod_fu_manager_initialized, &
   new_file_unit

!-----------------------------------------------------------------------
 
 implicit none

!-----------------------------------------------------------------------

! Module interface

 public :: &
   mod_landau_damping_test_constructor, &
   mod_landau_damping_test_destructor,  &
   mod_landau_damping_test_initialized, &
   test_name,        &
   test_description, &
   coeff_init

 private

!-----------------------------------------------------------------------

 ! public members
 character(len=*), parameter   :: test_name = "landau_damping"
 character(len=100), protected :: test_description(5)

 logical, protected :: mod_landau_damping_test_initialized = .false.

 ! private members
 real(wp), allocatable :: &
   k_wave(:),    & !< spatial wave number(s)
   vth(:),       & !< thermal velocity
   !> \f$ \left(\frac{v_{\rm th}^\alpha}{v_{\rm th}^1}\right)^2 \f$
   beta(:),      &
   amplitude(:)

 ! Input file ----------
 character(len=*), parameter :: &
   test_input_file_name = 'landau_damping.in'

 character(len=*), parameter :: &
   this_mod_name = 'mod_landau_damping_test'

 ! Input namelist, to be read from turbulent_channel.in
 namelist /input/ &
   k_wave, vth, amplitude

!-----------------------------------------------------------------------

contains

!-----------------------------------------------------------------------

 subroutine mod_landau_damping_test_constructor(d)
  integer, intent(in) :: d(2)

  integer :: fu, ierr
  character(len=*), parameter :: &
    this_sub_name = 'constructor'

   !Consistency checks ---------------------------
   if( (mod_messages_initialized.eqv..false.) .or. &
          (mod_kinds_initialized.eqv..false.) .or. &
     (mod_fu_manager_initialized.eqv..false.) ) then
     call error(this_sub_name,this_mod_name, &
                'Not all the required modules are initialized.')
   endif
   if(mod_landau_damping_test_initialized.eqv..true.) then
     call warning(this_sub_name,this_mod_name, &
                  'Module is already initialized.')
   endif
   !----------------------------------------------

   ! Allocate the module variables
   allocate( k_wave(d(1)) , amplitude(d(1)) , vth(d(2)) , beta(d(2)) )

   ! Read input file
   call new_file_unit(fu,ierr)
   open(fu,file=trim(test_input_file_name), &
      status='old',action='read',form='formatted',iostat=ierr)
    if(ierr.ne.0) call error(this_sub_name,this_mod_name, &
      'Problems opening the input file')
    read(fu,input)
   close(fu,iostat=ierr)

   beta = ( vth/vth(1) )**2

   write(test_description(1),'(a,i1,a,i1)') &
     'Landau_damping, space dimension ',d(1), &
                ', velocity dimension ',d(2)
   write(test_description(2),'(a,*(e12.6,:," , "))') &  
     '       space wave numbers: ',k_wave
   write(test_description(3),'(a,*(e12.6,:," , "))') &  
     '  perturbation amplitudes: ',amplitude
   write(test_description(4),'(a,*(e12.6,:," , "))') &  
     '       thermal velocities: ',vth
   write(test_description(5),'(a,*(e12.6,:," , "))') &  
     '      corresponding betas: ',beta

   mod_landau_damping_test_initialized = .true.
 end subroutine mod_landau_damping_test_constructor

!-----------------------------------------------------------------------
 
 subroutine mod_landau_damping_test_destructor()
  character(len=*), parameter :: &
    this_sub_name = 'destructor'
   
   !Consistency checks ---------------------------
   if(mod_landau_damping_test_initialized.eqv..false.) then
     call error(this_sub_name,this_mod_name, &
                'This module is not initialized.')
   endif
   !----------------------------------------------

   deallocate( k_wave , vth , beta , amplitude )

   mod_landau_damping_test_initialized = .false.
 end subroutine mod_landau_damping_test_destructor

!-----------------------------------------------------------------------

 !> Multivariate Gaussian profile with perturbations
 !!
 !! Denoting the space and velocity dimensions by \f$d_x\f$ and
 !! \f$d_v\f$, respectively, we have
 !! \f{displaymath}{
 !!  f(\underline{x},\underline{v}) = \left[\prod_{\alpha=1}^{d_v}
 !!  \sqrt{2\pi}\, v_{\rm th}^\alpha\right]^{-1}
 !!  \left( 1 - \sum_{\alpha = 1}^{d_x}
 !!  a^{\alpha}\sin(k^{\alpha}x^{\alpha}) \right)
 !!  \exp\left(-\frac{1}{2}\sum_{\alpha = 1}^{d_v}
 !!  \left(\frac{v^\alpha}{v_{\rm th}^\alpha} \right)^2 \right)
 !! \f}
 pure function coeff_init(x,v) result(f0)
  real(wp), intent(in) :: x(:,:), v(:,:)
  real(wp) :: f0(size(x,2))

  real(wp), parameter :: pi = 3.1415926535897932384626433832795028_wp
  integer :: l
  real(wp) :: ncoeff

  ncoeff = 1.0_wp/product(sqrt(2.0_wp*pi)*vth)
  do l=1,size(f0)
    f0(l) = ( 1.0_wp - sum(amplitude*sin(k_wave*x(:,l))) ) &
           * exp( -0.5_wp*sum((v(:,l)/vth)**2) )
  enddo
  f0 = ncoeff * f0 ! normalization

 end function coeff_init

!-----------------------------------------------------------------------

end module mod_landau_damping_test

